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Abstract 

We present a highly efficient molecular dynamics scheme for calculating the con- 
centration profile of dopants implanted in group-IV alloy, and III-V zinc blende 
structure materials. Our program incorporates methods for reducing computational 
overhead, plus a rare event algorithm to give statistical accuracy over several orders 
of magnitude change in the dopant concentration. The code uses a molecular dy- 
namics (MD) model, instead of the binary collision approximation (BCA) used in 
implant simulators such as TRIM and Marlowe, to describe ion-target interactions. 
Atomic interactions are described by a combination of 'many-body' and screened 
Coulomb potentials. Inelastic energy loss is accounted for using a Firsov model, 
and electronic stopping is described by a Brandt-Kitagawa model which contains 
the single adjustable parameter for the entire scheme. Thus, the program is easily 
extensible to new ion-target combinations with the minimum of tuning, and is pre- 
dictive over a wide range of implant energies and angles. The scheme is especially 
suited for calculating profiles due to low energy, large angle implants, and for situ- 
ations where a predictive capability is required with the minimum of experimental 
validation. We give examples of using our code to calculate concentration profiles 
and 2D 'point response' profiles of dopants in crystalline silicon, silicon-germanium 
blends, and gallium-arsenide. We can predict the experimental profile over five or- 
ders of magnitude for (100) and (110) channeling and for non-channeling implants 
at energies up to hundreds of keV. 

Key words: Molecular Dynamics, Ion Implant, Binary Collision, SIMS, 
Semiconductors . 



1 Introduction 

The continuing quest for greater processor performance demands ever smaller 
device sizes. The effort to realize these ultra shallow junction devices has re- 
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suited in current industry trends, such as the use of low energy, high mass, 
high dose, and large angle implants, to create abrupt dopant profiles. The 
experimental measurement of such profiles is challenging, as effects that are 
negligible at high implant energies become increasingly important as the im- 
plant energy is lowered. For example, the measurement of dopant profiles by 
secondary ion mass spectrometry (SIMS) is problematic for very low energy 
(less than 10 keV) implants, due to limited depth resolution of measured 
profiles. Also, refining SIMS protocols to obtain profiles for new ion-target 
combinations, e.g. In, Sb, or N implants or Sii^Ge^ targets is not a trivial 
problem. 

The use of computer simulation as an alternative method to determine dopant 
profiles is well established. Binary collision approximation (BCA) codes have 
traditionally been used, however such simulations become unreliable at low 
ion energies[l,2]. The BCA approach breaks down when multiple collisions 
or collisions between moving atoms become significant, or when the crystal 
binding energy is of the same order as the energy of the ion. Such problems 
are clearly evident when one attempts to use the BCA to simulate channeling 
in semiconductors; here the interactions between the ion and target are neither 
binary nor collisional in nature, rather they occur as many simultaneous soft 
interactions which steer the ion down the channel. 

A more accurate, alternative to the BCA, is the use of molecular dynamics 
(MD) simulation to calculate ion trajectories[3,4]. However, the computational 
cost of traditional MD simulations precludes the calculation of the thousands 
of ion trajectories necessary to construct a dopant profile. Here we present 
a highly efficient MD-based scheme, that is optimized to calculate the con- 
centration profiles of ions implanted into semiconductors. The algorithms are 
incorporated into our implant modeling molecular dynamics code [5], REED- 
MDQ. Our program has previously been demonstrated to describe the low 
dose (zero damage) implant of As, B, and P ions with energies in the sub 
MeV range into crystalline Si in (100), (110), and non-channeling directions, 
and also into amorphous Si [4-6]. We have now extended our model to any 
ion species, and to other diamond crystal substrates, such as C, Ge, SiC, 
Sii-zGe-j;, and GaAs. A model for ion induced damage has also been added to 
the program so that high dose implants can be simulated. 



2 Molecular Dynamics Model 

The basis of the molecular dynamics model is a collection of empirical poten- 
tial functions that describe interactions between atoms and give rise to forces 

1 Named for 'Rare Event Enhanced Domain following Molecular Dynamics'. 
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between them. In addition to the classical interactions described by the poten- 
tial functions, the interaction of the ion with the target electrons is required 
for ion implant simulations, as this is the principle way in which the ion loses 
energy. Another necessary ingredient is a description of the target material 
structure, including thermal vibration amplitudes. 

Interactions between target atoms are modeled by derivatives of the many- 
body potential developed by Tersoff[7,8]. ZBL 'pair specific' screened Coulomb 
potentials[9,10] are used to model interactions for common ion-target combi- 
nations. For other combinations, the ZBL 'universal' potential is used. The 
'universal' potential is also used to describe the close-range repulsive part of 
the Tersoff potentials. 

We include energy loss due to inelastic collisions, and energy loss due to elec- 
tronic stopping as two distinct mechanisms. The Firsov model [11] is used to 
describe the loss of kinetic energy from the ion due to momentum transfer 
between the electrons of the ion and target atom. We implement this using a 
velocity dependent pair potential, as derived by Kishinevskii[12]. 

A modified Brandt-Kitagawa[13] model, that involves both global and local 
contributions to the electronic stopping is used for electronic energy loss[4,6]. 
This model contains the single fitted parameter in our scheme, r°, the 'aver- 
age' one electron radius of the target material experienced by the ion. This 
is adjusted to account for oscillations in the Z\ dependence of the electronic 
stopping cross-section. The parameter is fit once for each ion-target combina- 
tion and is then valid for all ion energies and incident directions. By using a 
realistic stopping model, with the minimum of fitted parameters, we obtain a 
greater transferability to the modeling of implants outside the fitting set. 

In the calculations presented here, the target is a {100} diamond crystal with a 
surface oxide layer. The oxide structure was obtained from annealing a periodic 
Si0 2 sample with the density constrained to that estimated for grown surface 
oxide[14]. Thermal vibrations of atoms are modeled by displacing atoms from 
their lattice sites using a Debye model. For high dose implants, the accumu- 
lation of damage within the target is described by a simple Kinchin- Pease [15] 
model. Target properties are either species dependent, e.g., local electron den- 
sity, or are obtained by interpolation from known values for single element 
materials, e.g., lattice constant and Debye temperature. 



3 Efficient Molecular Dynamics Algorithms 

We apply a combination of methods to increase the efficiency of this specific 
type of simulation. Infrequently updated neighbor lists[16,17] are employed 
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to minimize the time spent in force calculations. The paths of the atoms 
are integrated using Verlet's algorithm [17], with a variable timestep that is 
dependent upon both kinetic and potential energy of atoms [16]. 

It is infeasible to calculate dopant profiles by full MD simulation, as com- 
putational requirements scale approximately as u 4 , where u is the initial ion 
velocity. We have developed a modified MD scheme which is capable of pro- 
ducing accurate dopant profiles with a much smaller computational overhead. 
We continually create and destroy target atoms, to follow the domain of the 
substrate that contains the ion. Material is built in front of the ion, and de- 
stroyed in its wake. Hence, the ion experiences the equivalent of a complete 
crystal, but the cost of the algorithm is only 0(u). 

To further improve efficiency, we use three other approximations. The mov- 
ing atom approximation [3] is used to reduce the number of force calculations. 
Atoms are divided into two sets; those that are 'on' have their positions inte- 
grated, and those that are 'off' are stationary. At the start of the simulation, 
only the ion is turned on. Some of the 'off' atoms will be used in the force 
calculations and will have forces assigned to them. If the resultant force ex- 
ceeds a certain threshold, the atom is turned on. We use two thresholds in 
our simulation; all atoms interacting directly with the ion are turned on im- 
mediately (zero threshold), and other atoms are turned on if the force exceeds 
8.0X10" 9 N. 

For high ion velocities, we do not need to use a many-body potential to main- 
tain the stable diamond lattice; a pair potential is sufficient, as only repulsive 
interactions are significant. Hence, at a certain ion velocity we switch from 
the complete many-body potential to a pair potential approximation for the 
target atom interactions. We make a further approximation for still higher 
ion energies, where only ion-target interactions are significant in determining 
the ion path. This approximation, termed the recoil interaction approxima- 
tion's], brings the MD scheme close to many BCA implementations. The 
major difference between the two approaches is that the ion path is obtained 
by integration, rather than by the calculation of asymptotes, and that multiple 
interactions are, by the nature of the method, handled in the correct manner. 



4 Rare Event Algorithm 

A typical dopant profile in a crystalline semiconductor consists of a near- 
surface peak followed by an almost exponential decay over several orders of 
magnitude in concentration. If we attempt to directly calculate a statistically 
significant dopant concentration at all depths of the profile we will have to 
run many ions that are stopped near the peak for every one ion that stops in 
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the tail, and most of the computational effort will not enhance the accuracy 
of the profile. 

In order to remove this redundancy, we employ an 'atom splitting' scheme[19] 
to increase the sampling in the deep component of the concentration profile. 
At certain splitting depths in the material, each ion is replaced by two ions, 
each with a statistical weighting of half that prior to splitting. Each split ion 
trajectory is run separately, and the weighting of the ion is recorded along with 
its final depth. As the split ions experience different environments (material is 
built in front of the ion, with random thermal displacements), the trajectories 
rapidly diverge from one another. Due to this scheme, we can maintain the 
same number of virtual ions moving at any depth, but their statistical weights 
decrease with depth. During a typical simulation, 1,000 implanted ions are 
split to yield around 10,000 virtual ions. 
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Fig. 1. The estimated uncertainty in the calculated dopant profile due to 2 keV As 

lxlO 12 cm" 2 (7,0) implant into Si with 1 unit cell of surface oxide, for the same 

number of initial ions (1,000), (a) with and (b) without rare event enhancement. 

We estimate the uncertainty in the calculated dopant profiles by dividing the 
final ion depths into 10 sets. A depth profile is calculated from each set using a 
histogram of 100 bins, and the standard deviation of the distribution of the 10 
concentrations for each bin is taken as the uncertainty. Fig. 1 demonstrates the 
effectiveness of the scheme, by comparing profiles obtained with and without 
atom splitting over five orders of magnitude. We estimate that the rare event 
algorithm reduces CPU time by a factor of 90 when calculating profiles over 
3 orders of magnitude, and by a factor of 900 when calculating a profile over 
5 orders of magnitude. 



5 Results and Discussion 

First, we give 2D profiles for low dose, low energy profiles to show scattering 
and surface effects. We then give examples of ID profiles produced by sim- 
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ulations, and compare to SIMS data[20,21]. All simulations were run with a 
target temperature of 300 K, and a beam divergence of 1.0° was assumed. 
Each profile was constructed from 1,000 ions, with the splitting depths up- 
dated every 25 ions, and a domain of 3x3x3 unit cells was used. The direction 
of the incident ion beam is specified by the angle of tilt, 9°, from normal and 
the azimuthal angle 0°, as (6,<f>). In the case of the low energy (< 10 keV) 
implants, we assume one unit cell thickness of surface oxide; for other cases 
we assume three unit cells of oxide at the surface. For the low energy implants, 
we have calculated profiles over a change of five orders of magnitude in con- 
centration; for the higher energy implants we calculate profiles over 3 orders 
of magnitude. A dose of lxlO 12 cm -2 (zero damage) is used unless otherwise 
noted. 
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2. Calculated 2D dopant profile due to 0.2 keV B (0,0) implant into Si. 
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Fig. 3. Calculated 2D dopant profiles due to 2 keV As (0,0) implants into Si, with 
(a) 1 and (b) 3 unit cells of surface oxide. 

2D profiles are shown projected onto the plane normal to the surface and con- 
taining the zero degree azimuth. This makes it easy to differentiate between 
major channeling directions; the (100) channel is vertical, and the four (110) 
channels appear at angles of 35° from vertical. Fig. 2 shows the result of an 
ultra-low energy implant into Si. Although the implant is in the (100) direc- 
tion, this channel is closed at such low ion energies and the only channeling 
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Fig. 4. Calculated dopant profiles for 10 keV B and As (0,0) implants into Sii_ x Ge x . 

occurs in the (110) direction. This demonstrates the need to have a 'universal' 
electronic stopping model, rather than a model tuned for a particular channel- 
ing direction. The effect of the amount of surface disorder is shown in Fig. 3. 
Increasing the thickness of the surface layer leads to more ions being scattered 
into the larger (110) channel, and hence gives a far deeper tail to the profile. 
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Fig. 5. Calculated and experimental[20] dopant profiles due to Zn 3xl0 13 cm~ 2 
implants into GaAs. 
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Fig. 6. The calculated and experimental [21] dopant profiles due to 150 keV Si 
3xl0 13 cm -2 implants into GaAs. 
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Fig. 7. The calculated and experimental[21] dopant profiles due to 300 keV Sc 
3xl0 13 cm~ 2 implants into GaAs. 

There is increasing interest in the use of SiGe as a replacement for Si cur- 
rently used in CMOS technology, due to its higher switching speed[22]. Fig. 4 
shows the effect of Ge concentration on profiles from B and As implants into 
Sii-^Ge-j; targets. The trend is clearly for shallower profiles with increasing Ge 
concentration, but this is extremely non-linear; the difference between x = 
and x = 0.2 profiles is greater than the difference between x = 0.2 and x = 0.8 
profiles. The remaining figures show the calculated concentration profiles of 
several ion species implanted under various conditions into GaAs substrates, 
and comparison with available SIMS data. The results of the REED calcu- 
lations show good agreement with the experimental data, demonstrating the 
accuracy of our model and its transferability to many ion-target combinations 
and implant conditions. 



6 Conclusions 



In summary, we have developed a restricted MD code to simulate the ion 
implant process and calculate 'as implanted' dopant profiles. This gives us the 
accuracy obtained by time integrating atom paths, with an efficiency far in 
excess of full MD simulation. 

The scheme described here gives a viable alternative to the BCA approach. 
Although it is still more expensive computationally, it is sufficiently fast to 
be used on modern desktop computer workstations. The method has two ma- 
jor advantages over the BCA approach: (i) Our MD model consists only of 
standard empirical potentials developed for bulk semiconductors and for ion- 
solid interactions. The only fitting is in the electronic stopping model, and 
this involves only one parameter per ion-target combination. We believe that 
by using physically based models for all aspects of the program, with the 
minimum of fitting parameters, we obtain good transferability to the mod- 
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eling of implants outside of our fitting set. (ii) The method does not break 
down at the low ion energies necessary for production of the next generation 
of semiconductor technology; it gives the correct description of multiple, soft 
interactions that occur both in low energy implants, and higher energy chan- 
neling. Hence our method remains truly predictive at these low ion energies, 
whilst the accuracy of the BCA is in doubt. 



Acknowledgements 

This work was performed under the auspices of the United States Department 
of Energy. 



References 



[I] M. T. Robinson and I. M. Torrens, Phys. Rev. B 12, 5008 (1974). 

[2] W. Eckstein, Computer Simulation of Ion-Solid Interactions (Springer- Verlag, 
New York, 1991) pp30-32. 

[3] D. E. Harrison, Jr., in Critical Reviews in Solid State and Materials Sciences, 
edited by J. E. Greene (CRC, Boca Raton, 1988), Vol. 14, Suppl. 1. 

[4] D. Cai, N. Gr0nbech-Jensen, C. M. Snell, and K. M. Beardmore, Phys. Rev. B 
54, 17147 (1996). 

[5] K. M. Beardmore and N. Gr0nbech- Jensen, Phys. Rev. E 57, 7278 (1998); 
jhttpj / /bifrost. lanl.gov/~reed/reed.shtm] 

[6] D. Cai, CM. Snell, K.M. Beardmore, and N. Gr0nbech-Jensen, Int. J. Mod. 
Phys. C 9, 459 (1998). 

[7] J. Tersoff, Phys. Rev. B 38, 9902 (1988); J. Tersoff, Phys. Rev. B 39, 5566 
(1989). 

[8] R. Smith, Nucl. Instr. and Meth. B 67, 335 (1992); M. Sayed, Nucl. Instr. and 
Meth. B 102, 218 (1995). 

[9] J. F. Ziegler, J. P. Biersack, and U. Littmark, The Stopping and Range of Ions 
in Solids (Pergamon Press, New York, 1985). 

[10] Potentials calculated from the ZBL tabulated electron densities, using a code 
written by E. Morvan. 

[II] O. B. Firsov, Zh. Eksp. Teor Fiz. 36, 1517 (1959) [Sov. Phys. JETP 9, 1076 
(1959)]. 



9 



[12] L. M. Kishinevskii, Izv. Acad. Nauk. SSSR, Ser. Fiz. 26, 1410 (1962) [Bull. 
Acad. Sci. USSR, Phys Ser. 26, 1433 (1962)]; V. A. Elteckov, D. S. Karpuzov, 
Yu. V. Martynenko, and V. E. Yurasova, in Atomic Collision Phenomena in 
Solids, Ed. D. W. Palmer, M. W. Thompson, and P. D. Townsend, (North 
Holland, Amsterdam, 1970) p. 657. 

[13] W. Brandt and M. Kitagawa, Phys. Rev. B 25, 5631 (1982). 

[14] A. Demkov, (private communication). 

[15] G. H. Kinchin and R. S. Pease, Rep. Prog. Phys. 18, 1 (1955). 

[16] K. M. Beardmore, Ph. D. Thesis, Loughborough University of Technology, 1995. 

[17] L. Verlet, Phys. Rev. 159, 98 (1967). 

[18] K. Nordlund, Comp. Mat. Sci. 3, 448 (1995). 

[19] J. M. Hammersley and D. C. Handscomb, Monte Carlo Methods (Methuen, 
London, 1964) p8. 

[20] R. G. Wilson, J. Appl. Phys. 53, 5641 (1982). 

[21] R. G. Wilson and V. R. Deline, Appl. Phys. Lett. 37, 793 (1980). 

[22] W. C. Holton, Solid State Tech. 40, 119 (November 1997). 



10 



